MODULE md1_hermite_expansion

   USE md1_global_consts
   USE md1_global_types
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: md1_build_Ecoefs

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE md1_build_Ecoefs(batch,xyz,Iang,Jang,Etij)

      IMPLICIT NONE
      TYPE(md1_prim_batch), INTENT(IN)  :: batch(:)
      INTEGER(INTK),        INTENT(IN)  :: xyz, Iang, Jang
      REAL(REALK),          INTENT(OUT) :: Etij(:,:,:,:)

      CALL md1_E_two_term_recurrence(batch,xyz,SIZE(batch),Iang,Jang,Etij)

!      REAL(REALK), ALLOCATABLE :: Etij_maw(:,:,:,:)
!      INTEGER(INTK) :: a,b,c,d


!INTEGER(INTK) :: i,j,t
!
!      ! Nakajima
!      IF (Iang >= Jang) THEN
!         CALL md1_build_Ecoef1(batch,xyz,SIZE(batch),Iang,Jang,Etij)
!      ELSE
!         CALL md1_build_Ecoef2(batch,xyz,SIZE(batch),Iang,Jang,Etij)
!      END IF
!
!return

      ! Watson
!      a = SIZE(Etij,1)
!      b = SIZE(Etij,2)
!      c = SIZE(Etij,3)
!      d = SIZE(Etij,4)
!      ALLOCATE (Etij_maw(a,b,c,d))
!      CALL md1_E_two_term_recurrence(batch,xyz,SIZE(batch),Iang,Jang,Etij_maw)

!do t = 1, b
!do i = 1, c
!do j = 1, d
!   IF (t > i+j-1) CYCLE
!   print *, t-1,i-1,j-1, Etij(:,t,i,j), " === ", Etij_maw(:,t,i,j)
!end do
!end do
!end do
!      IF (SUM(Etij) /= SUM(Etij_maw)) STOP 'md1_build_Ecoefs'

   END SUBROUTINE md1_build_Ecoefs

!-------------------------------------------------------------------------------

   SUBROUTINE md1_E_two_term_recurrence(batch,xyz,np,Iang,Jang,Etij)

!     o Generate E-coefficients (Hermite expansion coefficients) 
!       using the two-term recurence relations

      IMPLICIT NONE
      TYPE(md1_prim_batch), INTENT(IN)  :: batch(:)
      INTEGER(INTK),        INTENT(IN)  :: xyz, np, Iang, Jang
      REAL(REALK),          INTENT(OUT) :: Etij(:,0:,0:,0:)

      INTEGER(INTK) :: i,j,t, k
      REAL(REALK) :: p2(np)
      REAL(REALK) :: PA(np)
      REAL(REALK) :: PB(np)
      REAL(REALK) :: tmp(np)

      p2(:) = half / batch(1:np)%exp
      PA(:) = batch(1:np)%PA(xyz)
      PB(:) = batch(1:np)%PB(xyz)

      Etij(:,0:,0:,0:) = zero

      ! t   = 0
      ! i+j = 0
      IF (xyz == 3) THEN
         Etij(:,0,0,0) = batch(:)%KAB
      ELSE
         Etij(:,0,0,0) = one
      END IF

      DO k = 0, Iang+Jang
         ! t = 0
         ! i+j = k
         DO i = 0, Iang
            DO j = 0, Jang
               IF ( i+j /= k ) CYCLE
               IF (i>0) Etij(:,0,i,j) = PA(:)*Etij(:,0,i-1,j) + Etij(:,1,i-1,j)
               IF (j>0) Etij(:,0,i,j) = PB(:)*Etij(:,0,i,j-1) + Etij(:,1,i,j-1)
            END DO
         END DO
         ! t > 0
         ! i+j = t+k
         DO t = 1, Iang+Jang
            DO i = 0, Iang
               DO j = 0, Jang
                  IF ( i+j /= t+k ) CYCLE
                  tmp(:) = zero
                  IF ( i > 0 ) tmp(:) = tmp(:) + DBLE(i) * Etij(:,t-1,i-1,j)
                  IF ( j > 0 ) tmp(:) = tmp(:) + DBLE(j) * Etij(:,t-1,i,j-1)
                  Etij(:,t,i,j) = one / DBLE(t) * p2(:) * tmp(:)
               END DO
            END DO
         END DO
      END DO

   END SUBROUTINE md1_E_two_term_recurrence

!-------------------------------------------------------------------------------

   SUBROUTINE md1_build_Ecoef1(batch,xyz,NPrim,Iang,Jang,Etij)

      IMPLICIT NONE
      TYPE(md1_prim_batch), INTENT(IN)  :: batch(:)
      INTEGER(INTK),        INTENT(IN)  :: xyz, NPrim, Iang, Jang
      REAL(REALK),          INTENT(OUT) :: Etij(:,0:,0:,0:)

      INTEGER(INTK) :: I, J, t, JTemp
      REAL(REALK) :: ExpHalf(NPrim)
      REAL(REALK) :: PA(NPrim)
      REAL(REALK) :: PB(NPrim)
print *, iang, jang

!     o Generate E-coefficients (Hermite expansion coefficients) 
!       with the three-term recurence relation

      ExpHalf(:) = half / batch(1:Nprim)%exp
      PA(:) = batch(1:Nprim)%PA(xyz)
      PB(:) = batch(1:Nprim)%PB(xyz)

      IF (xyz == 3) THEN
         Etij(:,0,0,0) = batch(:)%KAB
      ELSE
         Etij(:,0,0,0) = One 
      END IF

      IF (Iang == 0) RETURN
      DO I = 1, Iang
         JTemp = I
         IF (I == Iang) JTemp = Jang
         DO J = 0, JTemp

            ! t = 0
            t = 0
            IF (I+J <= 1) THEN
               Etij(:,t,i,j) = PA(:) * Etij(:,t,i-1,j)
            ELSE
               Etij(:,t,i,j) = PA(:) * Etij(:,t,i-1,j) + Etij(:,t+1,i-1,j)
            END IF

            ! t >= 1
            DO t = 1, I+J
               IF (t == I+J) THEN
                  Etij(:,t,i,j) = ExpHalf(:) * Etij(:,t-1,I-1,J)
               ELSE IF (t == I+J-1) THEN
                  Etij(:,t,i,j) = ExpHalf(:) * Etij(:,t-1,I-1,J)   &
                                        + PA(:) * Etij(:,t,i-1,j)
               ELSE
                  Etij(:,t,i,j) = ExpHalf(:) * Etij(:,t-1,I-1,J)  &
                          + PA(:) * Etij(:,t,i-1,j)              &
                          + DBLE(t+1)      * Etij(:,t+1,i-1,j)
               END IF
            END DO

            IF (I /= j) THEN

               ! t = 0
               t = 0
               IF (I+J <= 1) THEN
print *, ubound(Etij,1)
print *, ubound(Etij,2)
print *, ubound(Etij,3)
print *, ubound(Etij,4)
print *, ubound(pb)
print *, t,i,j
                  Etij(:,t,j,i) = PB(:) * Etij(:,t,J,I-1)
               ELSE
                  Etij(:,t,j,i) = PB(:) * Etij(:,t,J,I-1) + Etij(:,t+1,J,I-1)
               END IF

               ! t >= 1
               DO t = 1, I+J
                  IF (t == I+J) THEN
                     Etij(:,t,j,i) = ExpHalf(:) * Etij(:,t-1,J,I-1)
                  ELSE IF (t == I+J-1) THEN
                     Etij(:,t,j,i) = ExpHalf(:) * Etij(:,t-1,J,I-1) &
                               + PB(:) * Etij(:,t,J,I-1)
                  ELSE
                     Etij(:,t,j,i) = ExpHalf(:) * Etij(:,t-1,J,I-1) &
                                     + PB(:) * Etij(:,t,J,I-1)  &
                                       + DBLE(t+1)      * Etij(:,t+1,J,I-1)
                     END IF
               END DO

            END IF

         END DO
      END DO

   END SUBROUTINE md1_build_Ecoef1

!-------------------------------------------------------------------------------

   SUBROUTINE md1_build_Ecoef2(batch,xyz,NPrim,Iang,Jang,Etij)

      IMPLICIT NONE
      TYPE(md1_prim_batch), INTENT(IN)  :: batch(:)
      INTEGER(INTK),        INTENT(IN)  :: xyz, NPrim, Iang, Jang
      REAL(REALK),          INTENT(OUT) :: Etij(:,0:,0:,0:)

      INTEGER(INTK) :: I, J, t, JTemp
      REAL(REALK) :: ExpHalf(NPrim)
      REAL(REALK) :: PA(NPrim)
      REAL(REALK) :: PB(NPrim)

!     o Generate E-coefficients (Hermite expansion coefficients) 
!       with the three-term recurence relation

      ExpHalf(:) = half / batch(1:Nprim)%exp
      PA(:) = batch(1:Nprim)%PA(xyz)
      PB(:) = batch(1:Nprim)%PB(xyz)

      IF (xyz == 3) THEN
         Etij(:,0,0,0) = batch(:)%KAB
      ELSE
         Etij(:,0,0,0) = One 
      END IF

      IF (Jang == 0) RETURN
      DO J = 1, Jang
         JTemp = J
         IF (J == Jang) JTemp = Iang
         DO I = 0, JTemp

            ! t = 0
            t = 0
            IF (I+J <= 1) THEN
               Etij(:,t,i,j) = PB(:) * Etij(:,t,I,j-1)
            ELSE
               Etij(:,t,i,j) = PB(:) * Etij(:,t,I,j-1) + Etij(:,t+1,I,j-1)
            END IF

            ! t >= 1
            DO t = 1, I+J
               IF (t == I+J  ) THEN
                  Etij(:,t,i,j) = ExpHalf(:) * Etij(:,t-1,I,j-1)
               ELSE IF (t == I+J-1) THEN
                  Etij(:,t,i,j) = ExpHalf(:) * Etij(:,t-1,I,j-1)  &
                              + PB(:)    * Etij(:,t,I,j-1)
               ELSE
                  Etij(:,t,i,j) = ExpHalf(:) * Etij(:,t-1,I,j-1) &
                                 + PB(:)    * Etij(:,t,I,j-1) &
                                  + DBLE(t+1)      * Etij(:,t+1,I,j-1)
               END IF
            END DO

            IF (I /= j) THEN

               ! t = 0
               t = 0
               IF (I+J <= 1) THEN
                  Etij(:,t,j,i) = PA(:) * Etij(:,t,J-1,I)
               ELSE
                  Etij(:,t,j,i) = PA(:) * Etij(:,t,J-1,I) + Etij(:,t+1,J-1,I)
               END IF

               ! t >= 1
               DO t = 1, I+J
                  IF (t == I+J  ) THEN
                     Etij(:,t,j,i) = ExpHalf(:) * Etij(:,t-1,J-1,I)
                  ELSE IF (t == I+J-1) THEN
                     Etij(:,t,j,i) = ExpHalf(:) * Etij(:,t-1,J-1,I) &
                                  + PA(:)    * Etij(:,t,J-1,I)
                  ELSE
                     Etij(:,t,j,i) = ExpHalf(:) * Etij(:,t-1,J-1,I) &
                                  + PA(:)    * Etij(:,t,J-1,I) &
                                   + DBLE(t+1)      * Etij(:,t+1,J-1,I)
                  END IF
               END DO

            END IF

         END DO
      END DO

   END SUBROUTINE md1_build_Ecoef2

!-------------------------------------------------------------------------------

END MODULE md1_hermite_expansion

